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Abstract 



Processing the signáis received on an array 
of sensors for the location of the emitter is of great 
enough interest to have been treated under many 
special case assumptions. The general problem 
considers sensors with arbitrary locations and 
arbitrary directional characteristics (gain phase 
polarization) in a noise/interference environment of 
arbitrary covariance matrix. This report is 
concerned first with the múltiple emitter aspect of 
this problem and second with the generality of 
solution. A description is given of the múltiple 
signal classification (MUSIC) algorithm, which 
provides asymptotically unbiased estimates of 1) 
number of incident wavefronts present; 2) directions 
of arrival (DOA) (or emitter locations); 3) strengths 
and cross correlations among the incident 
waveforms; 4) noise/interference strength. Examples 
and comparisons with methods based on máximum 
likelihood (ML) and máximum entropy (ME), as well 
as conventional beamforming are. included. An 
example of its use as a múltiple frequency estimator 
operating on time series is included. 

INTRODUCTION 

The Term Múltiple signal classification 
(MUSIC) is used to describe experimental and 
theoretical techniques involved in determining the 
parameters of múltiple wavefronts arriving at an 
antenna array from measurements made on the signáis 
received at the array elements. 

The general problem considers antennas with 
arbitrary locations and arbitrary directional 
characteristics (gain/phase/ polarization) in a 
noise/interference environment of arbitrary 
covariance matrix. The múltiple signal classification 
approach is described; it can be implemented as an 
algorithm to provide asymptotically unbiased 
estimates of 

1) number of signáis; 

2) directions of arrival (DOA); 

3) strengths and cross correlations among the 
directional waveforms; 



4) polarizations; 

5) strength of noise/interference. 

These techniques are very general and of 
wide application. Special cases of MUSIC are 

1) conventional interferometry; 

2) monopulse direction finding (DF), i.e., using 
múltiple colocated antennas; 

3) múltiple frequency estimation 

THE DATA MODEL 

The waveforms received at the M array 
elements are linear combinations of the D incident 
wavefronts and noise. Thus, the múltiple signal 
classification approach begins with the following 
model for characterizing the received M vector X as 
in 



=[a(e 1 )a(e 2 ).a(e D )] 













w 2 




+ 











Or 

X=AF+ W. (1) 

The incident signáis are represented in 
amplitude and phase at some arbitrary reference point 
(for instance the origin of the coordínate system) by 

the complex quantities F 1} F 2 , , F D . The noise, 

whether "sensed" along with the signáis or generated 
internal to the instrumentation, appears as the 
complex vector W. 

The elements of X and A are also complex in 
general. The a¡¡ are known functions of the signal 
arrival angles and the array element locations. That 
is, a¡j depends on the ¡th array element, its position 
relative to the origin of the coordínate system, and its 
response to a signal incident from the direction of the 
jth signal. The jth column of A is a 'mode' vector 
a( 6j) of responses to the direction of arrival 9j of the 
jth signal. 

In geometrical language, the measured X 
vector can be visualized as a vector in M dimensional 
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space. The directional mode vectors a( 6J) - a¡j for i = 
1, 2, ... M, i.e.. the columns of A, can also be se 
visualized. Equation (1) states that X is a particular 
linear combination of the mode vectors, the elements 
of Fare the coefficients of the combination. Note that 
the X vector is confined to the range space of A. That 
is, if A has two columns, the range space is no more 
than a two-dimensional subspace within the M space 
and X necessarily lies in the subspace Also note that 
a6, the continuum of all possible mode vectors, lies 
within the M space but is quite nonlinear. For help in 
visualizing this, see Fig. 1. For example, in an 
azimuth-only direction finding system, 8 will consist 
of a single parameter. In an azimuth/elevation/ range 
system, 9 will be replaced by 9, (|), r for example. In 
any case, a(9) is a vector continuum such as a "snake" 
(azimuth only) or a "sheet" (Az/El) twisting and 
winding through the M space. (In practice, the 
procedure by which the a(9) continuum is measured 
or otherwise established corresponds to calibrating 
the array. 



TME SIGNAL SU&SPACE 
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Fig. ! . Gcometric portrayal for tfirce-antenna case. 



In these geometrical terms (see Fig. 1), the 
problem of solving for the directions of arrival of 
múltiple incidents continuum with the range space of 
A. The range space of A is course, obtained from the 
measured data. The means of obtaining the range 
space and, necessarily, its dimensionality (the number 



D is incident signáis) follows. 
The S Matrix 

The MxM covariance matrix of the X vector 

is 



S=XX* =AFF*A + WW 



S= APA*+ X S 

Under the basic assumption that the incident 
signáis and the noise are uncorrelated. Note that the 
incident waveforms represented by the elements of F 

may be uncorrelated (the DxX matrix P =FF* is 
diagonal) or may contain completely correlated pairs 
(P is singular). In general, P will be "Merely" positive 
definite which reflects the arbitrary degrees of pairs - 
wise correlations oceurring between the incident 
waveforms. 

When the number of incidents wavefronts D 
is less than the number of array elements M, then 
APA* is singular; it has a rank less than M. Therefore 

IAPA*l = IS-^S o l=0 

(3) 

This equation is only satisfied with X equal 
to one of the eigenvalues of S in the metric of So. 
But, for A full rank and P positive definite, APA* 
must be nonnegative definite. Therefore X can only 
be the minimum eigenvalue A, min . Therefore, any 
measured S = XX* matrix can be written 

S*APA* + A mi „S , Xmin>o 

where X min is the smallest solution to IS - XS \ 
= 0. Note the special case wherein the elements of the 
noise vector W are mean zero, variance a 2 , in which 

Case, ^minSo = cj 2 I 

CALCULATING A SOLUTION 

The rank of APA* is D and can be determined 
directly from the eigenvalues of S in the metric of So, 
That is, in the complete set of eigenvalues of S in the 
metric of So» ^ m in will not always be simple. In fact, 
it oceurs repeated N = M - D times. This is true 
because the eigenvalues of S and those of S - ^ m i n S = 
APA* differ by X min in all cases. Since the minimum 
eigenvalue of APA* is zero (being singular), X min must 
oceur repeated N times. Therefore, the number of 
incident signáis estimator is 

"D=M-"N 

where = the multiplicity of ^ min (S, S ) and 
X m { B (S, So) is read A, m i„ of S in the metric of So-" (In 
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practice, one can expect that the múltiple X min will 
occur in a cluster rather than all precisely equal. The 
1 \spread" on this cluster decreases as more data is 
processed.) 

The Signal and Noise Subspaces 

The M eigenvectors of S in the metric of S 
must satisfy' Se¡ = A,¡ S e, if i = 1,2, M Since S = APA* 
+ ^minS , we have APA*e¡ ¡ = (X-, - X min )S e i Clearly, 
for each of the X-¡ that is equal to X m - ln there are N — 
we must have APA*e¡ = or A*e¡ = 0. That is, the 
eigenvectors associated with X min (S, S ) are 
orthogonal to the space spanned by the columns of A; 
the incident signal mode vectors! 

Thus we may justifiably refer to the N 
dimensional subspace spanned by the N noise 
eigenvectors as the noise subspace and the D 
dimensional subspace spanned by the incident signal 
mode vectors as the signal subspace: they are 
disjoint. 

THE ALGORITHM 

We now have the means to solve for the 
incident signal mode vectors. If E N is defined to be 
the M x N matrix whose columns are the N noise 
eigenvectors, and the ordinary Euclidean distance 
(squared) from a \ector Y to the signal subspace is d 2 - 
Y*E N E N Y. we can plot Vd for points along the a( 0) 
continuum as a function of 9. That is, 

1 



MU 



(6) 



a*(e)E N E* N a(0) 



(However, the a ( 9) continuum may intersect the D 
dimensional signal subspace more than D times; 
another unresolvable situation occurring only for the 
case of múltiple incident signáis — a type II 
ambiguity.) It is clear from the expression that 
MUSIC is asymptotically unbiased even for múltiple 
incident wavefronts because S is asymptotically 
perfectly measured so that E N is also. a(6) does not 
depend on the data. 1 , Once the directions of arrival of 
the D incident signáis have been found, the A matrix 
becomes available and may be used to compute the 
parameters of the incident signáis. The solution for 
the P matrix is direct 1 and can be expressed in terms 
of (S - ^mi n S ) and A. That is, since APA* = S - 

P= (A*A)"' A* (S-X mm S ) A (A*A) _1 

INCLUDING POLARIZATION 

Consider a signal arriving from a specific 
direction 9 D . Assume that the array is not diverse in 



polarization; Le., all elements are identically 
polarized, say, vertically. Certainly the DF system 
will be most sensitive to vertically polarized energy, 
completely insensitive to horizontal and partially 
sensitive to arbitrarily polarized energy. The array is 
only sensitive to the vertically polarized component 
of the arriving energy. 

For a general or polarizationally diverse 
array, the mode vector corresponding to the direction 
O depends on the signal polarization. A vertically 
polarized signal will induce one mode vector and 
horizontal another, and right-hand circular (RHC) 
still another. 

Recall that signal' polarization can be 
completely characterized by a single complex number 
q. We can 1 'observe' how the mode vector changes as 
the polarization parameter q for the emitter changes 
at the specific direction 9 Q . It can be oroven that as q 
changes through all possible polarizations, the mode 
vector sweeps out a two-dimensional "polarization 
subspace." Thus, only two independent mode vectors 
spanning the polarization subspace for the direction 
O are needed to represent any emitter polarization q 
at direction O . The practical embodiment of this is 
that only the mode vectors of two emitter 
polarizations need be ejaculated or kept in store for 
direction O in order to solve for emitter polarizations 
where only one was needed to solve for DOA in a 
system with an array that was not polarizationally 
diverse. 

These arguments lead to an equation similar 
to (6) for P( 0) but including the effects of 
polarization diversity among the array elements. 

1 



Pmu(6) = 



/r 



X. 



I(e) 

:(e) 



E N E* N [a x (0)a y (0)] 



(8) 



where a x ( 0) and a y ( 0) are the two continua 
corresponding to, for example, separately taken x and 
y linear incident wave- 
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y y. y» y 
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a*(e)S" 1 a(e) 

1 
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front polarizations. The eigenvector corresponding to 
^min in (8) provides the polarization parameter q since 
it is of the form [1 q] T . 

THE ALGORITHM 

In summary, the steps of the algorithm (see 
Fig. 2) are: 

StepO: collect data form S; 

Step 1: calcúlate eigenstructure of S in metric of 

s ; 

Step 2: decide number of signáis D; (5); 

Step 3: evalúate P M u(&) versus (9) (6) or (8): 

Step 4: pick D peaks of Pmu( 0)1 

Step 5: calcúlate remaining parameters; (7). 

The above steps have been implemented in 
several forms to verify and evalúate the principies 
and basic performance. Field tests have been 
conducted using actual receivers, arrays, and múltiple 
transmitters. The results of these tests have 
demonstrated the potential of this approach for 
handling múltiple signáis in practical situations. 
Performance results are being prepared for 
presentation in another paper. 

COMPARISON WITH OTHER METHODS 

In comparing MUSIC with ordinary 
beamforrning (BF), máximum likelihood (ML), and 
máximum entropy (ME), the following expressions 
were used. See Figs. 3 and. 4. 

P BF (0) = a*(9)Sa(0) 
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Fig. 3. Example of azimuth-only DF performance. 
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where c is a column of S . The beamformer 
expression calculates for plotting the power one 
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would measure at the output of a beamformer 
(summing the array element signáis after inserting 
delays appropriate to steer or look in a specific 
direction) as a function of the direction. 

Pml(8) calculates the log likelihood function under 
the assumptions that X' is a mean zero, multivariate 
Gaussian and that there is only a single incident 
wavefront present. For 

múltiple incident wavefronts, P M ( 9) becomes 

which implies a D dimensional search (and 

plot!). 

P ME (8) is based on selecting one of the M 
array elements as a "reference" and attempting to 
find weights to be applied to the remaining M - 1 
received signáis to permit their sum with a MMSE fit 
to the reference. Since there are M possible 
references, there are M generally different Pme( 0) 
obtained from the M possible column selections from 
S' 1 . In the comparison plots, a particular reference 
was consistently selected. 

An example of the completely general 
MUSIC algorithm applied to a problem of steering a 
múltiple feed parabolic dish antenna is shown in Fig. 
5. sin x/x pencil beamshapes skewed slightly off 
boresight are assumed for the element patterns. Since 
the six antennas are essentially colocated, the DF 
capacity arises out of the antenna beam pattern 
diversity. The computer was used to simúlate the 
"noisy" S matrix that would arise in practice for the 
conditions desired and then to subject it to the 
MUSIC algorithm. Fig. 5 shows how three directional 
signáis are distinguished and their polarizations 
estimated even though two of the arriving signáis are 
highly similar (90 percent correlated). 

The application of MUSIC to the estimation 
of the frequencies of múltiple sinusoids (arbitrary 
amplitudes and phases) for a very limited duration 
data sample is shown in Fig. 6. The figure suggests 
that, even though there was no actual noise included, 
the rounding of the data samples to six decimal digits 
has already destroyed a significant portion of the 
information present in the data needed to resolve the 
several frequencies) 

SUMMARY AND CONCLUSION 

As this paper was being prepared, the works 
of Gething [1] and Davies [2] were discovered, 



offering a part of the solution discussed here but in 
terms of simultaneous equations and special linear 
relationships without recourse to eigenstructure. 
However, the geometric significance of a vector space 
setting and the interpretation of the S matrix 
eigenstructure was missed. More recent work by 
Reddi [3] is also along the lines of the work presented 
here though limited to uniform, collinear arrays of 
omnidirectional elements and also without clear 
utilization of the entire noise subspace. Ziegenbein 
[4] applied the same basic concept to time series 
spectral analysis referring to it as a Karhunen-Loeve 
transform though treating aspects of it as "ad hoc." 
El-Behery and MacPhie [5] and Capón [6] treat the 
uniform collinear array of omnidirectional elements 
using the máximum likelihood method. Pisarenko [7] 
also treats time series and addresses only the case of 
a full complement of sinusoids; i.e., a one- 
dimensional noise subspace. 

The approach presented here for múltiple signal 
classification is very general and of wide application. 
The method is interpretable in terms of the geometry 
of complex M spaces wherein the eigenstructure of 
the measured S matrix plays the central role. MUSIC 
provides asymptotically unbiased estimates of a 
general set of signal parameters, approaching Cramer- 
Rao accuracy bound, MUSIC models the data as tr.e 
sum of point source emissions and noise rather than 
vr.t convolution of an all pole transfer function 
driver, by a white noise (i.e., autoregressive 
modeling, máximum entropy; or maximizing a 
probability under the assumption that the A" vector is 
zero mean, Gaussian (máximum likelihood for 
Gaussian data), In geometric terms MUSIC minimizes 
the distance from the a(9) continuum to the signal 
subspace whereas máximum likelihood minimizes a 
weighted combination all component distances. 

No assumptions have been made about array 
geometry. The array elements may be arranged in a 
regular or irregular pattern and may differ or be 
identical in directional characteristics 
(amplitude/phase) provided their polarization 
characteristics are all identical. The extensión to 
include general polarizationally diverse antenna 
arrays will be more completely described in a 
sepárate paper. 
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